INFLUENCE OF NATURAL CONVECTION AND THERMAL RADIATION 
MULTI-COMPONENT TRANSPORT IN MOCVD REACTORS 
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The influence of Grashof and Reynolds number in Metal Organic Chemical Vapor (MOCVD) reactors is being 
investigated under a combined empirical/numerical study. As part of that research, the deposition of Indium 
Phosphide in an MOCVD reactor is modeled using the computational code CFD-ACE. The model includes the 
effects of convection, conduction, and radiation as well as multi-component diffusion and multi-step surface/gas 
phase chemistry. The results of the prediction are compared with experimental data for a commercial reactor and 
analyzed with respect to the model accuracy. 

Keywords: Chemical Vapor Deposition (CVD), MOCVD, Computational Fluid Dynamics (CFD), Indium Phosphide, 
Radiation, Richardson Number 


5 * y-3 


-y / 


^ <y 


1. INTRODUCTION 


Thermal Predictions 


Metalorganic Chemical Vapor Deposition (MOCVD) is a 
common technique for the growth of III- V and II- VI 
compound semiconductors and alloys. The uniformity and 
quality of these materials is tightly coupled to the composition 
and temperature distribution at the growth interface. As such, 
the understanding and control of heat and mass transport in 
MOCVD reactors is critical to the production of high quality 
materials. 

The mechanisms of heat and mass transport during MOCVD 
are currently being investigated under a joint project between 
CFD Research Corporation and NASA Langley Research 
Center (LaRC). The research is a combined numerical/ 
experimental study to investigate the influence of radiation and 
natural convection in MOCVD systems over a range of 

Richardson (Gr/Re 2 ) numbers. 

2. NUMERICAL MODEL 

The computational fluid dynamics code CFD-ACE is being 
used in the current study to model and analyze the chemistry 
and heat and mass transfer for selected MOCVD reactors. 
The code is a three-dimensional control volume Navier-Stokes 
code with surface and gas phase chemistry. Several advanced 
features essential for modeling MOCVD systems have been 
incorporated into the code. These include models for non-gray 
radiation, Soret effects, multi-component diffusion, and 
advanced surface chemistry models using CHEMKIN. 

3. RESULTS TO DATE 

The numerical studies performed thus far can be place in three 
categories: 1) thermal validation, 2) convection analysis, and 
3) deposition predictions. These are discussed in more detail 
below. 


Thermal studies have been performed using data obtained by 
LaRC for pure gas flow in an MOCVD sled reactor. The 
apparatus used for these experiment is located in the Chemical 
Vapor Deposition Facility for Reactor Characterization at 

NASA Langley Research Center (Figure l) 1 . The reactor has a 
circular inlet section that feeds the reactants into a rectangular 
duct in which is mounted a fused silica sled containing a 
graphite susceptor. The graphite susceptor is heated by an 
external Radio Frequency (RF) induction coil (not shown in 
Figure 1), wound around the outside of the silica test section. 
The reactor was operated at various flow conditions using 
pure nitrogen, helium, and hydrogen (Table 1). During the 
experiment, the temperature along the sides and top wall of 

the reaction chamber was measured using an infrared camera^. 
These data were subsequently used for model validation* . 
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Figure 1 . Experimental MOCVD Reactor at LaRC. 
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Table 1 . Pure Gas Thermal/Flow Experiments 
Conducted at LaRC 


Gr/Re 2 = 27 

Hydrogen 

Helium 

Nitrogen 

Flow Rate 

9.9 Ipm 

N/A 

9.9 1pm 

Gr 

200 

N/A 

9690 

Re 

2.7 

N/A 

18.9 


Gr/Re 2 = 41.5 

Hydrogen 

Helium 

Nitrogen 

Flow Rate 

8 1pm 

8 1pm 

8 Ipm 

Gr 

200.0 

171.0 

9690 

Re 

2.2 

2.0 

15.3 


Gr/Re 2 = 166 

— 

Hydrogen 

Helium 

Nitrogen 

Flow Rate 

4 1pm 

4 Ipm 

4 1pm 

Gr 

200 

171.0 

9690 

Re 

1.1 

1.0 

7.6 


Gr/Re 2 = 664 

Hydrogen 

Helium 

Nitrogen 

Flow Rate 

N/A 

N/A 

2 1pm 

Gr 

N/A 

N/A 

9690 

Re 

N/A 

N/A 

3.8 


Gr/Re 2 = 
2655 

Hydrogen 

Helium 

Nitrogen 

Flow Rate 

N/A 

N/A 

1 1pm 

Gr 

N/A 

N/A 

9690 

Re 

N/A 

N/A 

1.9 


All the test cases in Table l were re-run using the 
numerical code. The model domain for these 
simulations consists of the entire reactor as shown in 
Figure 1, i.e., the chamber, the graphite susceptor, 


the fused silica walls, and the fused silica susceptor 
holder. The gray radiation model assumption was 
used. The temperature on the reactor wall was 
determined as part of the solution, rather than 
specified a priori. The induction heating was 
modeled by fixing the graphite to 873K, as set in 
the experiment. The results of the simulation 
compare well with the thermal data 3 . For the range 
of Richardson numbers simulated, radiation and 
convection were both found to play significant roles 
in the heat transport in the reactor. Figure 2 shows a 
comparison of the empirical data with the 
numerical predictions for Nitrogen at Ri = 41.5 
where Ri is the Richardson number, defined as Gr/ 

Re . (The “dips” in the experimental date are due to 
the RF coil obstructing the IR camera view). The 
agreement between the experimental data and the 
numerical model (e.g. to within 10 degrees for the 
wall temperature above the substrate) is considered 
to be satisfactory for first order simulations of the 
reactor environment. However, for very accurate 
analyses, additional physics, such as the non-gray 
radiation and the transient effects of surface 
emissivities need to be included. These mechanisms 
have been shown to be significant in commercial 
reactors and, as such, need to be precisely 
understood and modeled. 

Silica Wail Temperatures For Nitrogen (Rj=41.5) 
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Figure 2. CFD-ACE: Thermal Prediction Valida- 
tion for Nitrogen Flow in a CVD Reactor. 
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Flow Simulations 


Reaction Simulation 


The parametric set of experimental test cases 
performed at LaRC were also analyzed numerically 
in terms of recirculation flows as a function of 
Richardson number. These parametric cases 
reconfirm the dependence of recirculation on the 
Richardson number and the existence of both 
longitudinal and transverse rolls for this reactor 

design 5 ' 9 . Figure 3 show a comparison of the 
relative location and magnitude of the three 
dimensional recirculation zones inside the reactors 
as indicated by an iso-surface of reverse flow equal 
to one cm/sec. The results indicate the presence of 
recirculation. It was determined numerically that 
changing the reactor tilt could reduce the 
recirculation zones, however, since the simulations 
are for a non-reacting system, it is not possible to 
draw a correlation between the magnitude of the 
recirculation and the effect it would have on 
uniformity of deposition. This leads to the next 
phase of the study in which a reacting system is 
modeled. 



Figure 3. CFD-ACE: Recirculations Patterns for 
Pure Gases in an MOCVD as a Function 
of Ri #. 


The deposition of Indium Phosphide (InP) from the 
precursors Phosphine (PH 3 ), Trimethylindium 
(In(CH 3 ) 3 ), and Monomethylindium (InCH 3 ) was 
selected as the initial material for modeling. This 
system was selected because of the importance of 
InP as a semiconductor material and because of the 
availability of experimental data from the University 

of Virginia 10 . Furthermore, the data for this study 
were obtained using a commercial Crystal 
Specialties Model 425 MOCVD reactor (Figure 4) 
very similar in design to the test apparatus at LaRC . 



Figure 4. MOCVD Reactor at the University of 
Virginia. 

The reaction set assumed for this system and the 
corresponding rate equations are shown in Tables 2 
& 3 respectively 10 : 

Table 2. Reaction Rates 


Surface 


In(CH 3 )^ 

+ PH^ —> InP(s) + 3CH 4 

(la) 

InCH-j + PH^ -> lnP(s) + CH 4 + H 2 
GasPhase 

In(CH^)^ — » InCH^ + 2CH^ 

(lb) 


Table 3. Reaction Equations 


^InP “ ' 

f -E./RT -E,/RT] 

| A j (TM I]e +A 3 [MMI]e j 

(2a) 

W T M1 = ' 

-E./RT -E 5 /RT ] 

-A j [TMlje -A 2 [TMI]e 

(2b) 
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f -E-,/RT -E./RT) 

“MMI = A 2 [TMI]e ‘ -A 3 [MMI]e (2c) 

f -E./RT -E./RT) 

cbpH^ = j- A i[ TMI ]e -A 3 [MMI]e l (2d) 

f -E./RT -E,/RT) 

d) CH4 = j3Aj[TMI]e +A 3 (MMI]e \ (2e) 

f -E ? /RT } 

d) CH = j2A 2 [TMI]e ‘ (2f) 

f -E-./RT) 

= WMMIJe | (2g) 

The activations energies for the reaction rates are 
based on the experimental measurements of 
Buchan, et al 11 , and the pre-exponential factors 
were determined based on calibration with the 
empirical results 10 . 


Thermal Conductivity (W / m - K): 

k = 9.6336 x 10' 3 + 3.453 x 10' 4 • T 

- 1.413e-08*T 2 (5c) 

The run conditions used for the simulation are the 
same as UVA InP experimental run # 83 as listed 
below in Table 7. 

Table 7. Run Conditions 

Total Flow: 

7720 seem (cm 3 /min) (1.1479 e-5 kg/s)* 

TMI (Trimethylindium) inlet concentration: 

467e-6 mole fraction (4.844e-4 mass fraction) ** 
TMI flow rate: 0.048 seem 
PH3 inlet concentration: 

3.88e-3 mole fraction (0.06179 mass fraction) ** 
PH3 flow rate: 30 seem 


Table 4. Activation Energies 


A, = 5 x 10 5 

1/sec 

(3a) 

14 

A 2 = 1 x 10 

1/sec 

(3b) 

A, = lx I0 9 

1 /sec 

(3c) 


Table 5. Pre-exponential Coefficients 


= 7.6 x l() 7 J/(kgmol) 

(4a) 

= 1 .7 x 1 (}**.)/( kg mol ) 

(4b) 

= 1 .51 x l()^J/(kginol ) 

(4c) 


With the exception of the diffusivities, which are 
computed for each species based on kinetic theory, 
the physical properties of the mixture were assumed 
to be those of the predominate hydrogen. These 
properties were computed as follows 10 : 

Table 6. Properties 

Viscosity (kg / m - sec): 

H = 2.907 x 10‘ 6 +2. 173 x 10" 8 • T (5a) 

- 4.9167 x 10 I2 *T 2 

Specific Heat (J/kg - K): 

Cp =1.491 x 10 4 - 1.644 *T+ 1.709 x 10' 3 *T 2 (5b) 


Substrate temperature: 600°C (873 °K) 

Pressure: 760 torr ( 1 .01 3e4 nt/m 3 or 1 .0 atm) 

* Based on a density of 0.0893 kg/m 3 at T = 

273°K and pressure = 101300 nt/m 2 (1 atm) 

** Based on a mixture molecular weight of 2. 1 349 

The Reynolds Number (Re = ULp/ji) for this 
simulation is 18. This is based on a length scale of 3 
centimeters which is the channel height at the 
leading edge of the susceptor. The viscosity is 
8.983e-6 kg/m-s, based on the inlet temperature 
(300K) and the viscosity formula in Table 6 The 
corresponding mixture density, at 1.0 atmosphere 
and 300 K, is 0.0862 kg/m 3 , as computed using the 
ideal gas law. The reference velocity, U, is .06255 
m/s based on the above density, the cross-sectional 
area at the susceptor leading edge (2.743e-2 m 2 ), 
and the mass flow rate (1.1479e-5 kg/s). 

The Grashof Number (Gr = g (3 [T hoI - T cold ] L 3 /v 2 ) 
for the simulation is 16.000. The T hot is the 
substrate temperature of 873 °K and the T cold is the 
inlet temperature of 300 °K. For an ideal gas, the 
thermal coefficient of expansion ((3) is equal to 1/T. 
For this case, |3 is approximated as a constant, (3=1/ 
T ref . where T re f is set to the susceptor temperature 
(873/K). The length scale, L, and the kinematic 
viscosity v = ju/p (1.042 e-3 m 2 /s) are the same as 
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used for the Reynolds number and g is one earth’s 
gravity. 

The corresponding Richardson number (Gr/Re 2 ) is 
49 

Simulation Results 

Flow Pattern and Temperature Distribution: 

Figure 5 shows the predicted flow pattern over the 
sled. At the Richardson number of 49, the flow 
pattern is dominated by forced convection, with 
minor recirculation zones above the leading edge of 
the substrate. 


CFD-ACE: InP Case 83 Ri = 49 
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Figure 5. CFD-ACE Centerline Flow Pattern 
(Velocity Vectors) for UVA Case 83. 

The corresponding centerline temperature contours 
are as shown in Figure 6. 


CFD-ACE: InP Case 83 Ri = 49 
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Figure 6. CFD-ACE 

Centerline 

Temperature 

Contours. 




Predicted Deposition Rate: A full 3-D cut-away 

view (similar to the angle in Figure 1) of the model 
prediction is shown in Figure 7. The deposition rate 
of InP on top of the sled is indicated by the grey 
flooded surface contours and the x-y plot in the 
lower right-hand of the figure. The measured 


deposition rate is included in the x-y plot as a dotted 
line. The deposition rate along the centerline 
exhibits a double peak near the leading edge of the 

substrate, as reported by Black 10 . This is caused by 
the dual surface deposition mechanisms ( 1 a and 1 b) 
listed in Table 2. Two-dimensional studies with a 
refined grid do a better job of reproducing the 
sharpness of these twin peaks. However, even grid 
refinement does not increase the average deposition 
rate predicted over the majority of substrate. In 
effect, the model is under predicting the deposition 
rate by nearly a factor of 2. 

CVD Growth of inP 
1 CFD-ACE Model 

| 


i 



||M- Case 33 ; > : : 

|jFU* : 49 1.0 Alffr 7720 seem 

Figure 7. CFD-ACE Prediction Deposition Rates 
and Selected Surface Temperatures for 
InP at 7720 seem at 1 atm (Ri = 49). 

In an attempt to identify the cause of the 
discrepancy between the measured and predicted 
growth rates, several parametric studies were 
conducted to identify the sensitivity of the growth 
rate to the model assumptions. These include grid 
refinement, inlet swirl, modification of the pre- 
exponential coefficients (by 10%), gas radiative 
absorption, and substrate temperature (by 30 K). 
None of these were shown to increase the deposition 
rate significantly. This indicates that the predicted 
growth is diffusion limited. As shown by Figure 8, 
the concentration of MMI next to the substrate is 
severely depleted, thus limiting the growth rate. 
This was confirmed by increasing the reaction rate 

in Eq. 3c from 10 9 to 10 20 with only a minor 
increase in the average growth rate. 
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CFD-ACE: InP Case 83 Ri = 49 



lnCH 3 Mass Fraction (Along Centerline} 
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Figure 8. CFD-ACE Centerline Mass Fraction 
Contours. 

The question therefore remains as to why the 
predicted growth rate of InP is less than the 
experimental rate. A very plausible explanation is 

provided by Black 10 who points out that the InP 
growth in the experiments is polycrystalline and that 
in the initial stage of growth, it is relatively sparse, 
as shown in Figure 9. This is in contrast with the 
model which computes growth rates in terms of kg/ 

m -s and assumes a perfect surface. (Black’s own 
simulation show a similar discrepancy between the 
experiments and predictions). 



Figure 9. SEM micrograph of InP deposition on 
fused-silica after a growth time of 2 
hours. Growth conditions the same as 
Run 32 (0. 1 atm, 4420 seem) 10 . 

The above discrepancies point to the need for both 
better models and well controlled experiments. 
Towards that end, the current model will be 
extended to include more rigorous chemistry models 
with more complex gas phase chemistry and surface 

reactions including surface adsorbed species 12 ' 14 


The material to be modeled with this advanced 
model will be the growth of GaAs, which is well 
characterized and for which there are high quality 
data available for comparison 15 ' 16 . 

4. SUMMARY/CONCLUSIONS 

A full 3-D CFD model of a commercial MOCVD 
reactor has been developed and applied to model 
the growth of Indium Phosphide. The thermal 
predictive capability of the model has previously 
been validated using empirical data for pure gas 
flow. The model predictions for InP deposition, 
based on a simplified set of reaction equations are 
within a factor of two of the measured deposition. 
The empirical data itself is shown to be highly 
dependent on the polycrystalline structure and 
varies experimentally by a factor of two over a 
growth period of five hours. Parametric numerical 
studies indicates that the predicted growth rate is 
diffusion limited and relatively insensitive to grid 
refinement, increases is substrate temperature, inlet 
swirl, gas absorptivity, and variations in the pre- 
experimental factors. 

The discrepancy between the predicted rate and the 
data is attributed to the polycrystalline nature of the 
growth surface and to the simplified chemistry 
mechanisms in the model. 

The next phase of the project will study the 
deposition process with more complex (and 
realistic) chemistry models that account for surface 
adsorbed species. Once these improvements have 
been sufficiently validated, the model will then be 
used to better understand the combined role of 
radiation/convection on the MOCVD deposition 
process. 
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